module TopolParser

using Printf
using Utils

export to_itp, to_rtp, get_atom_from_nr, get_nr_from_atom, gen_topol_from_file, resort!, append_top!
export gen_topol_atom_from_line, gen_topol_bond_from_line, gen_topol_angle_from_line, gen_topol_dihedral_from_line
export gen_topol_constraint_from_line, gen_topol_pair_from_line, gen_topol_exclustion_from_line
export TopolAtom, TopolBond, TopolConstraint, TopolPair, TopolAngle, TopolDihedral, TopolExclusion, Topol
export add_atom!, add_bond!, add_constraint!, add_pair!, add_angle!, add_dihedral!, add_exclusion!

struct TopolAtomtype
    name::String
    atnum_or_type::String
    mass::Float64
    charge::Float64
    ptype::String
    sigma::Float64
    epsilon::Float64
end

mutable struct TopolAtom
    nr::Int32
    type::String
    resnr::Int32
    resid::String
    atom::String
    cgnr::Int32
    charge::Float64
    mass::Float64       # -1: default mass
end

mutable struct TopolBond
    ai::TopolAtom
    aj::TopolAtom
    funct::Int32
    c0::Float64
    c1::Float64
end

mutable struct TopolPair
    ai::TopolAtom
    aj::TopolAtom
    funct::Int32
end

mutable struct TopolConstraint
    ai::TopolAtom
    aj::TopolAtom
    funct::Int32
    cs::Vector{Float64}
end

mutable struct TopolAngle
    ai::TopolAtom
    aj::TopolAtom
    ak::TopolAtom
    funct::Int32
    c0::Float64
    c1::Float64
end

mutable struct TopolDihedral
    ai::TopolAtom
    aj::TopolAtom
    ak::TopolAtom
    al::TopolAtom
    funct::Int32
    c0::Float64
    c1::Float64
    c2::Int32
end

mutable struct TopolExclusion
    atoms::Vector{TopolAtom}
end

mutable struct Topol
    atomtypes::Set{TopolAtomtype}
    moleculetype::String
    nrexcl::Int32
    atoms::Vector{TopolAtom}
    bonds::Vector{TopolBond}
    pairs::Vector{TopolPair}
    constraints::Vector{TopolConstraint}
    angles::Vector{TopolAngle}
    dihedrals::Vector{TopolDihedral}
    exclusions::Vector{TopolExclusion}
end

function get_atom_from_nr(atoms::Vector{TopolAtom}, nr)
    """从原子编号获取原子拓扑对象"""
    return filter(x -> x.nr == nr, atoms)[1]
end


function get_nr_from_atom(atoms, resnr, atname)
    """从原子拓扑对象获取原子编号"""
    return filter(x -> x.resnr == resnr && x.atom == atname, atoms)[1]
end

function add_atom!(atoms::Vector{TopolAtom}, atom_top; pos=-1)
    insert!(atoms, pos == -1 ? length(atoms) + 1 : pos, atom_top)
    resort!(atoms)
end

function add_bond!(bonds::Vector{TopolBond}, bond_top)
    push!(bonds, bond_top)
end

function add_constraint!(constraints::Vector{TopolConstraint}, constraint_top)
    push!(constraints, constraint_top)
end

function add_pair!(pairs::Vector{TopolPair}, pair_top)
    push!(pairs, pair_top)
end

function add_angle!(angles::Vector{TopolAngle}, angle_top)
    push!(angles, angle_top)
end

function add_dihedral!(dihedrals::Vector{TopolDihedral}, dihedral_top)
    push!(dihedrals, dihedral_top)
end

function add_exclusion!(exclusions::Vector{TopolExclusion}, exclusion_top)
    push!(exclusions, exclusion_top)
end

function gen_topol_from_file(topfile::AbstractString)::Topol
    println("Reading topology of $(basename(topfile))...")
    lines = strip.(readlines(topfile))
    filter!(x -> !occursin(r"^[;\#].*", x), lines)
    filter!(!isempty, lines)
    lines = (line -> split(line, ";")[1]).(lines)

    # topol items
    attypes = Set{TopolAtomtype}()
    mol = "MOL"
    nrexcl = 3
    atoms = Vector{TopolAtom}([])
    bonds = Vector{TopolBond}([])
    pairs = []
    constraints = []
    angles = []
    dihedrals = []
    exclusions = []

    cur_item = ""

    for line in lines
        if (m = match(r".*\[(.*)].*", line)) !== nothing
            cur_item = m.captures[1]
            cur_item = strip(cur_item)
        else
            if cur_item == "atomtypes"
                push!(attypes, gen_topol_attype_from_line(line))
            elseif cur_item == "moleculetype"
                mol = split(line)[1]
                nrexcl = parse(Int32, split(line)[2])
            elseif cur_item == "atoms"
                push!(atoms, gen_topol_atom_from_line(line))
            elseif cur_item == "bonds"
                push!(bonds, gen_topol_bond_from_line(atoms, line))
            elseif cur_item == "pairs"
                push!(pairs, gen_topol_pair_from_line(atoms, line))
            elseif cur_item == "constraints"
                push!(constraints, gen_topol_constraint_from_line(atoms, line))
            elseif cur_item == "angles"
                push!(angles, gen_topol_angle_from_line(atoms, line))
            elseif cur_item == "dihedrals"
                push!(dihedrals, gen_topol_dihedral_from_line(atoms, line))
            elseif cur_item == "exclusions"
                push!(exclusions, gen_topol_exclustion_from_line(atoms, line))
            else
                continue
            end
        end
    end
    resort!(atoms)
    finish_hint("Finished reading topology of " * mol * ".\n")
    return Topol(attypes, mol, nrexcl, atoms, bonds, pairs, constraints, angles, dihedrals, exclusions)
end

function gen_topol_attype_from_line(line::AbstractString)::TopolAtomtype
    paras = split(line)
    if length(paras) == 7
        name, atnum_or_type, mass, charge, ptype, sigma, epsilon = paras
    elseif length(paras) == 6
        name, mass, charge, ptype, sigma, epsilon = paras
        atnum_or_type = ""
    end
    mass = parse(Float64, mass)
    charge = parse(Float64, charge)
    sigma = parse(Float64, sigma)
    epsilon = parse(Float64, epsilon)
    return TopolAtomtype(name, atnum_or_type, mass, charge, ptype, sigma, epsilon)
end

function gen_topol_atom_from_line(line::AbstractString)::TopolAtom
    paras = split(line)
    if length(paras) == 8
        nr, type, resnr, resid, atom, cgnr, charge, mass = paras
        mass = parse(Float64, mass)
    else
        nr, type, resnr, resid, atom, cgnr, charge = paras
        mass = -1
    end
    nr = parse(Int32, nr)
    resnr = parse(Int32, resnr)
    cgnr = parse(Int32, cgnr)
    charge = parse(Float64, charge)
    return TopolAtom(nr, type, resnr, resid, atom, cgnr, charge, mass)
end

function gen_topol_bond_from_line(atoms::Vector{TopolAtom}, line::AbstractString)::TopolBond
    paras = split(line)
    ai, aj, funct = parse.(Int32, paras[1:3])
    ai = get_atom_from_nr(atoms, ai)
    aj = get_atom_from_nr(atoms, aj)
    if length(paras) > 3
        c0, c1 = parse.(Float64, paras[4:end])
    else
        c0 = -1
        c1 = -1
    end
    return TopolBond(ai, aj, funct, c0, c1)
end

function gen_topol_constraint_from_line(atoms::Vector{TopolAtom}, line::AbstractString)::TopolConstraint
    paras = split(line)
    ai, aj, funct = parse.(Int32, paras[1:3])
    ai = get_atom_from_nr(atoms, ai)
    aj = get_atom_from_nr(atoms, aj)
    if length(paras) > 3
        cs = parse.(Float64, paras[4:end])
    else
        cs = []
    end
    return TopolConstraint(ai, aj, funct, cs)
end

function gen_topol_pair_from_line(atoms::Vector{TopolAtom}, line::AbstractString)::TopolPair
    paras = split(line)
    ai, aj, funct = parse.(Int32, paras)
    ai = get_atom_from_nr(atoms, ai)
    aj = get_atom_from_nr(atoms, aj)
    return TopolPair(ai, aj, funct)
end

function gen_topol_angle_from_line(atoms::Vector{TopolAtom}, line::AbstractString)::TopolAngle
    paras = split(line)
    ai, aj, ak, funct = parse.(Int32, paras[1:4])
    ai = get_atom_from_nr(atoms, ai)
    aj = get_atom_from_nr(atoms, aj)
    ak = get_atom_from_nr(atoms, ak)
    if length(paras) > 4
        c0, c1 = parse.(Float64, paras[5:end])
    else
        c0 = -1
        c1 = -1
    end
    return TopolAngle(ai, aj, ak, funct, c0, c1)
end

function gen_topol_dihedral_from_line(atoms::Vector{TopolAtom}, line::AbstractString)::TopolDihedral
    paras = split(line)
    ai, aj, ak, al, funct = parse.(Int32, paras[1:5])
    ai = get_atom_from_nr(atoms, ai)
    aj = get_atom_from_nr(atoms, aj)
    ak = get_atom_from_nr(atoms, ak)
    al = get_atom_from_nr(atoms, al)
    if length(paras) > 5
        if funct == 1 || funct == 4 || funct == 9
            c0, c1, c2 = parse.(Float64, paras[6:end])
        elseif funct == 2
            c0, c1 = parse.(Float64, paras[6:end])
            c2 = -1
        else
            throw("Error: Dihedral parameters too much: $line")
        end
    else
        c0 = -1
        c1 = -1
        c2 = -1
    end
    return TopolDihedral(ai, aj, ak, al, funct, c0, c1, c2)
end

function gen_topol_exclustion_from_line(atoms::Vector{TopolAtom}, line::AbstractString)::TopolExclusion
    paras = split(line)
    atnums = parse.(Int32, paras)
    atoms = (x -> get_atom_from_nr(atoms, x)).(atnums)
    return TopolExclusion(atoms)
end

function shift_atoms!(top::Topol, offset::Int32)
    for atom in top.atoms
        if atom.nr > 0
            atom.nr += offset
            atom.cgnr += offset
        end
    end
    println("Atom numbers of ", top.moleculetype, " shifted by ", offset)
end

function shift_residues!(top::Topol, offset)
    # shift 残基编号
    for atom in top.atoms
        atom.resnr += offset
    end
    println("Residue numbers of ", top.moleculetype, " shifted by ", offset)
end

function Base.show(io::IO, top_atomtype::TopolAtomtype)
    @printf(io, "%-8s%6s%13.6f%12.6f%5s%18.6E%16.6E", top_atomtype.name, top_atomtype.atnum_or_type, 
    top_atomtype.mass, top_atomtype.charge, top_atomtype.ptype, top_atomtype.sigma, top_atomtype.epsilon)
end

function to_rtp(top_atom::TopolAtom)
    return @sprintf("%8s%6s%12.6f%5d", top_atom.atom, top_atom.type, top_atom.charge, top_atom.cgnr)
end

function Base.show(io::IO, top_atom::TopolAtom)
    # 12345678901234567890123456789012345678901234567890123456789
    #     1  HS14    1    FMC     HN2    1    0.245   1.0080
    @printf(io, "%7d%7s%7d%7s%8s%7d%12.6f", top_atom.nr, top_atom.type, top_atom.resnr, top_atom.resid, 
    top_atom.atom, top_atom.cgnr, top_atom.charge)
    if top_atom.mass != -1
        @printf(io, "%9.4f", top_atom.mass)
    end
end

function to_rtp(top_bond::TopolBond)
    s = @sprintf("%7s%7s", top_bond.ai.atom, top_bond.aj.atom)
    if top_bond.c0 != -1 && top_bond.c1 != -1
        s *= @sprintf("%13.4f%13.4e", top_bond.c0, top_bond.c1)
    end
    return s
end

function Base.show(io::IO, top_bond::TopolBond)
    # 12345678901234567890123456789012345678901234567890123456789
    #     1    2    2   0.0985   8.2455e+07
    @printf(io, "%7d%7d%9d", top_bond.ai.nr, top_bond.aj.nr, top_bond.funct)
    if top_bond.c0 != -1 && top_bond.c1 != -1
        @printf(io, "%13.4f%13.4e", top_bond.c0, top_bond.c1)
    end
end

function Base.show(io::IO, top_pair::TopolPair)
    # 12345678901234567890123456789012345678901234567890123456789
    #     1    5    1
    @printf(io, "%7d%7d%9d", top_pair.ai.nr, top_pair.aj.nr, top_pair.funct)
end

function Base.show(io::IO, top_constraint::TopolConstraint)
    # 12345678901234567890123456789012345678901234567890123456789
    #     1    5    1
    @printf(io, "%7d%7d%9d", top_constraint.ai.nr, top_constraint.aj.nr, top_constraint.funct)
    if top_constraint.cs != []
        (x -> @printf(io, "%13.4f", x)).(top_constraint.cs)
    end
end

function to_rtp(top_angle::TopolAngle)
    s = @sprintf("%7s%7s%7s", top_angle.ai.atom, top_angle.aj.atom, top_angle.ak.atom)
    if top_angle.c0 != -1 && top_angle.c1 != -1
        s *= @sprintf("%10.2f%9.2f", top_angle.c0, top_angle.c1)
    end
    return s
end

function Base.show(io::IO, top_angle::TopolAngle)
    # 12345678901234567890123456789012345678901234567890123456789
    #     1    5    1
    @printf(io, "%7d%7d%7d%9d", top_angle.ai.nr, top_angle.aj.nr, top_angle.ak.nr, top_angle.funct)
    if top_angle.c0 != -1 && top_angle.c1 != -1
        @printf(io, "%10.2f%9.2f", top_angle.c0, top_angle.c1)
    end
end

function to_rtp(top_dihedral::TopolDihedral)
    s = @sprintf("%7s%7s%7s%7s", top_dihedral.ai.atom, top_dihedral.aj.atom, top_dihedral.ak.atom, top_dihedral.al.atom)
    if top_dihedral.funct == 2
        if top_dihedral.c0 != -1 && top_dihedral.c1 != -1
            s *= @sprintf("%10.2f%9.2f", top_dihedral.c0, top_dihedral.c1)
        end
    elseif top_dihedral.funct == 4 || top_dihedral.funct == 1 || top_dihedral.funct == 9
        if top_dihedral.c0 != -1 && top_dihedral.c1 != -1 && top_dihedral.c2 != -1
            s *= @sprintf("%10.2f%9.2f%8d", top_dihedral.c0, top_dihedral.c1, top_dihedral.c2)
        end
    end
    return s
end

function Base.show(io::IO, top_dihedral::TopolDihedral)
    # 12345678901234567890123456789012345678901234567890123456789
    #     1    5    1
    @printf(io, "%7d%7d%7d%7d%9d", top_dihedral.ai.nr, top_dihedral.aj.nr, top_dihedral.ak.nr, 
    top_dihedral.al.nr, top_dihedral.funct)
    if top_dihedral.funct == 2
        if top_dihedral.c0 != -1 && top_dihedral.c1 != -1
            @printf(io, "%10.2f%9.2f", top_dihedral.c0, top_dihedral.c1)
        end
    elseif top_dihedral.funct == 4 || top_dihedral.funct == 1 || top_dihedral.funct == 9
        if top_dihedral.c0 != -1 && top_dihedral.c1 != -1 && top_dihedral.c2 != -1
            @printf(io, "%10.2f%9.2f%8d", top_dihedral.c0, top_dihedral.c1, top_dihedral.c2)
        end
    end
end

function to_rtp(top_exclusion::TopolExclusion)
    return join((x -> @sprintf("%7s", x.atom)).(top_exclusion.atoms))
end

function Base.show(io::IO, top_exclusion::TopolExclusion)
    # 12345678901234567890123456789012345678901234567890123456789
    #     1    5    1
    (x -> @printf(io, "%7d", x.nr)).(top_exclusion.atoms)
end

# TODO: 如果是末端残基, 就把Nxxx和Cxxx一起生成
function to_rtp(top, rtp_path, ff, exclude_N, exclude_C, atom_N, atom_C)
    # 首先把前后残基中的原子名加前缀
    for atom in top.atoms
        if atom in exclude_N && atom != atom_N
            if !startswith(atom.atom, "-")
                atom.atom = "-" * atom.atom
            end
        elseif atom in exclude_C && atom != atom_C
            if !startswith(atom.atom, "+")
                atom.atom = "+" * atom.atom
            end
        end
    end

    f = open(rtp_path, "w")
    write(f, "; rtp created by SuperMDA (https://gitee.com/supernova_bingbing/SuperMDA)\n")
    write(f, "; converted from top of $(top.moleculetype)\n\n")
    if ff == "amber"
        write(f, "[ bondedtypes ]\n")
        write(f, "; bonds  angles  dihedrals  impropers all_dihedrals nrexcl HH14 RemoveDih\n")
        write(f, "     1       1          9          4        1         3      1     0\n\n")
    elseif ff == "gromos"
        write(f, "[ bondedtypes ]\n; bonds  angles  dihedrals  impropers\n")
        write(f, "    2       2          1          2\n\n")
    else
        error_hint("Error: invalid forcefield, currently only support amber and gromos.\n")
        return
    end

    # 残基名
    write(f, "[ $(top.moleculetype) ]\n")

    # [ atoms ]字段：记录残基中每个原子的名称、类型和电荷、电荷组
    write(f, " [ atoms ]\n")
    for atom in top.atoms
        if !isempty(exclude_N) && atom in exclude_N
            write(f, "; $(to_rtp(atom))\t; previous residue\n")
        elseif !isempty(exclude_C) && atom in exclude_C
            write(f, "; $(to_rtp(atom))\t; next residue\n")
        else
            write(f, to_rtp(atom) * "\n")
        end
    end

    # [ bonds ]字段：原子间的连接信息
    write(f, " [ bonds ]\n")
    # amber力场保留-, gromos力场保留+
    filter!(bond -> !(all(atom -> atom in append!(exclude_N, exclude_C), [bond.ai, bond.aj])), top.bonds)
    if ff == "amber"
        filter!(bond -> !(bond.ai in exclude_C || bond.aj in exclude_C), top.bonds)
    elseif ff == "gromos"
        filter!(bond -> !(bond.ai in exclude_N || bond.aj in exclude_N), top.bonds)
    end
    for bond in top.bonds
        if !isempty(exclude_N) && (bond.ai in exclude_N && bond.aj in exclude_N)
            write(f, ";- " * to_rtp(bond) * "\n")
        elseif !isempty(exclude_C) && (bond.ai in exclude_C && bond.aj in exclude_C)
            write(f, ";+ " * to_rtp(bond) * "\n")
        else
            write(f, to_rtp(bond) * "\n")
        end
    end

    # [ angles ]字段：键角信息
    write(f, " [ angles ]\n")
    # 不要带超过1个多余原子的
    filter!(angle -> !(count(atom -> atom in append!(exclude_N, exclude_C), [angle.ai, angle.aj, angle.ak]) > 1), top.angles)
    for angle in top.angles
        if !isempty(exclude_N) && (angle.ai in exclude_N && angle.aj in exclude_N && angle.ak in exclude_N)
            write(f, ";- " * to_rtp(angle) * "\n")
        elseif !isempty(exclude_C) && (angle.ai in exclude_C && angle.aj in exclude_C && angle.ak in exclude_C)
            write(f, ";+ " * to_rtp(angle) * "\n")
        else
            write(f, to_rtp(angle) * "\n")
        end
    end

    # [ dihedrals ]字段：二面角信息
    write(f,  " [ dihedrals ]\n")
    # 9类型不要带多余原子的
    filter!(dihedral -> 
        !(count(atom -> atom in append!(exclude_N, exclude_C), [dihedral.ai, dihedral.aj, dihedral.ak, dihedral.al]) > 1), 
        top.dihedrals)
    for dihedral in top.dihedrals
        if dihedral.funct == 9 || dihedral.funct == 1
            if !isempty(exclude_N) && (dihedral.ai in exclude_N && dihedral.aj in exclude_N
                                    && dihedral.ak in exclude_N && dihedral.al in exclude_N)
                write(f,  ";- " * to_rtp(dihedral) * "\n")
            elseif !isempty(exclude_C) && (dihedral.ai in exclude_C && dihedral.aj in exclude_C
                                      && dihedral.ak in exclude_C && dihedral.al in exclude_C)
                write(f,  ";+ " * to_rtp(dihedral) * "\n")
            else
                write(f,  to_rtp(dihedral) * "\n")
            end
        end
    end

    # [ impropers ]字段：异常二面角信息
    write(f,  " [ impropers ]\n")
    for dihedral in top.dihedrals
        if dihedral.funct == 4 || dihedral.funct == 2
            if !isempty(exclude_N) && (dihedral.ai in exclude_N && dihedral.aj in exclude_N
                                    && dihedral.ak in exclude_N && dihedral.al in exclude_N)
                write(f,  ";- " * to_rtp(dihedral) * "\n")
            elseif !isempty(exclude_C) && (dihedral.ai in exclude_C && dihedral.aj in exclude_C
                                      && dihedral.ak in exclude_C && dihedral.al in exclude_C)
                write(f,  ";+ " * to_rtp(dihedral) * "\n")
            else
                write(f,  to_rtp(dihedral) * "\n")
            end
        end
    end

    # [ exclusions ]字段：定义哪些原子之间的非键相互作用忽略掉
    write(f,  " [ exclusions ]\n")
    # 超过一个的不要
    filter!(exclusion -> !(any((atom -> atom in exclude_N).(exclusion.atoms)) ||
                           any((atom -> atom in exclude_C).(exclusion.atoms))), top.exclusions)
    for exclusion in top.exclusions
        if !isempty(exclude_N) && all((x -> x in exclude_N).(exclusion))
            write(f,  ";- " * to_rtp(exclusion) * "\n")
        elseif !isempty(exclude_C) && all((x -> x in exclude_C).(exclusion))
            write(f,  ";+ " * to_rtp(exclusion) * "\n")
        else
            write(f, to_rtp(exclusion) * "\n")
        end
    end

    close(f)
end

function to_itp(top, file_name)
    f = open(file_name, "w")
    write(f, "; Created by SuperMDA (https://gitee.com/supernova_bingbing/SuperMDA)\n")
    # 输出原子类型
    if !isempty(top.atomtypes)
        write(f, "\n[ atomtypes ]\n; name   at.num      mass       charge   ptype     sigma (nm)    epsilon (kJ/mol)\n")
        (attype -> write(f, "$attype\n")).(top.atomtypes)
    end
    # 输出残基名
    write(f, "\n[ moleculetype ]\n; Name   nrexcl\n")
    write(f, "  $(top.moleculetype)    $(top.nrexcl)\n")
    # 输出原子编号
    write(f, "\n[ atoms ]\n;    nr   type  resnr  resid    atom   cgnr      charge     mass\n")
    (atom -> write(f, "$atom\n")).(top.atoms)
    # 输出键长编号
    write(f, "\n[ bonds ]\n;    ai     aj    funct           c0           c1\n")
    (bond -> write(f, "$bond\n")).(top.bonds)
    # 输出 pair 编号
    if !isempty(top.pairs)
        write(f, "\n[ pairs ]\n;    ai     aj    funct  ;  all 1-4 pairs but the ones excluded in GROMOS itp\n")
        (pair -> write(f, "$pair\n")).(top.pairs)
    end
    # 输出 constraint 编号
    if !isempty(top.constraints)
        write(f, "\n[ constraints ]\n")
        (constraint -> write(f, "$constraint\n")).(top.constraints)
    end
    # 输出键角编号
    write(f, "\n[ angles ]\n;    ai     aj     ak    funct     angle       fc\n")
    (angle -> write(f, "$angle\n")).(top.angles)
    # 输出 proper 二面角编号
    write(f, "\n[ dihedrals ]; proper\n;    ai     aj     ak     al    funct     phase       kd      pn   ; funct = 9 or 1\n")
    for dihedral in top.dihedrals
        if dihedral.funct == 1 || dihedral.funct == 9
            write(f, "$dihedral\n")
        end
    end
    # 输出 improper 二面角编号
    write(f, "\n[ dihedrals ]; improper\n;    ai     aj     ak     al    funct      ksi         k           ; funct = 2\n")
    write(f, ";    ai     aj     ak     al    funct     phase       kd      pn   ; funct = 4\n")
    for dihedral in top.dihedrals
        if dihedral.funct == 2 || dihedral.funct == 4
            write(f, "$dihedral\n")
        end
    end
    # 输出 exclusion 编号
    if !isempty(top.exclusions)
        write(f, "\n[ exclusions ]\n")
        (exclusion -> write(f, "$exclusion\n")).(top.exclusions)
    end
    close(f)
end

function resort!(atoms::Vector{TopolAtom})
    for (i, atom) in enumerate(atoms)
        atom.nr = i
        atom.cgnr = i
    end
    first_resnr = atoms[1].resnr
    (atom -> atom.resnr -= (first_resnr - 1)).(atoms)
end

function append_top!(top::Topol, another_top::Topol; overlap_file::String="", current_first::Bool=true)
    if overlap_file != ""
        lines = readlines(overlap_file)
        over_a1 = atrange2atlist(lines[1])
        over_a2 = atrange2atlist(lines[2])
        over_b1 = atrange2atlist(lines[3])
        over_b2 = atrange2atlist(lines[4])
        if length(over_a1) != length(over_a2) || length(over_b1) != length(over_b2)
            error_hint("Group B1 and B2 must be in one-to-one correspondence!!!\n")
            return
        end
    else
        over_a1 = []
        over_a2 = []
        over_b1 = []
        over_b2 = []
    end

    # Combine atom types
    println("Merging atom types of $(another_top.moleculetype) with $(top.moleculetype) ...")
    union!(top.atomtypes, another_top.atomtypes)

    # Combine bonds part
    println("Appending bonds of $(another_top.moleculetype) to $(top.moleculetype) ...")
    new_bonds_list = []
    # 处理主体
    for bond in top.bonds
        if length(over_b1) != 0 && length(over_b2) != 0
            if current_first                                # 排除主体B1-B1, 主体A1-B1替换成A1-B2
                if !(bond.ai.nr in over_b1) || !(bond.aj.nr in over_b1)
                    if bond.ai.nr in over_b1
                        bond.ai = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== bond.ai.nr)])
                    end
                    if bond.aj.nr in over_b1
                        bond.aj = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== bond.aj.nr)])
                    end
                    push!(new_bonds_list, bond)
                end
            else                                            # 排除主体A1-B1
                if !(bond.ai.nr in over_b1) && !(bond.aj.nr in over_b1)
                    push!(new_bonds_list, bond)
                end
            end
        else
            push!(new_bonds_list, bond)
        end
    end
    # 处理客体
    for bond in another_top.bonds
        if length(over_a1) != 0 && length(over_a2) != 0
            if current_first                                                # 排除客体A2-A2
                if !(bond.ai.nr in over_a2) && !(bond.aj.nr in over_a2)
                    push!(new_bonds_list, bond)
                end
            else                                                            # 排除客体A2-A2, 客体A2-B2替换为A1-B2
                if !(bond.ai.nr in over_a2) || !(bond.aj.nr in over_a2)
                    if bond.ai.nr in over_a2
                        bond.ai = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== bond.ai.nr)])
                    end
                    if bond.aj.nr in over_a2
                        bond.aj = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== bond.aj.nr)])
                    end
                    push!(new_bonds_list, bond)
                end
            end
        else
            push!(new_bonds_list, bond)
        end
    end
    top.bonds = new_bonds_list

    # Combine pairs part
    println("Appending pairs of $(another_top.moleculetype) to $(top.moleculetype) ...")
    new_pairs_list = []
    # 处理主体
    for pair in top.pairs
        if length(over_b1) != 0 && length(over_b2) != 0
            if current_first                                # 排除主体B1-B1, 主体A1-B1替换成A1-B2
                if !(pair.ai.nr in over_b1) || !(pair.aj.nr in over_b1)
                    if pair.ai.nr in over_b1
                        pair.ai = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== pair.ai.nr)])
                    end
                    if pair.aj.nr in over_b1
                        pair.aj = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== pair.aj.nr)])
                    end
                    push!(new_pairs_list, pair)
                end
            else                                            # 排除主体A1-B1
                if !(pair.ai.nr in over_b1) && !(pair.aj.nr in over_b1)
                    push!(new_pairs_list, pair)
                end
            end
        else
            push!(new_pairs_list, pair)
        end
    end
    # 处理客体
    for pair in another_top.pairs
        if length(over_a1) != 0 && length(over_a2) != 0
            if current_first                                                # 排除客体A2-A2
                if !(pair.ai.nr in over_a2) && !(pair.aj.nr in over_a2)
                    push!(new_pairs_list, pair)
                end
            else                                                            # 排除客体A2-A2, 客体A2-B2替换为A1-B2
                if !(pair.ai.nr in over_a2) || !(pair.aj.nr in over_a2)
                    if pair.ai.nr in over_a2
                        pair.ai = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== pair.ai.nr)])
                    end
                    if pair.aj.nr in over_a2
                        pair.aj = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== pair.aj.nr)])
                    end
                    push!(new_pairs_list, pair)
                end
            end
        else
            push!(new_pairs_list, pair)
        end
    end
    top.pairs = new_pairs_list

    # Combine constraints part
    println("Appending constraints of $(another_top.moleculetype) to $(top.moleculetype) ...")
    new_constraints_list = []
    # 处理主体
    for constraint in top.constraints
        if length(over_b1) != 0 && length(over_b2) != 0
            if current_first                                # 排除主体B1-B1, 主体A1-B1替换成A1-B2
                if !(constraint.ai.nr in over_b1) || !(constraint.aj.nr in over_b1)
                    if constraint.ai.nr in over_b1
                        constraint.ai = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== constraint.ai.nr)])
                    end
                    if constraint.aj.nr in over_b1
                        constraint.aj = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== constraint.aj.nr)])
                    end
                    push!(new_constraints_list, constraint)
                end
            else                                            # 排除主体A1-B1
                if !(constraint.ai.nr in over_b1) && !(constraint.aj.nr in over_b1)
                    push!(new_constraints_list, constraint)
                end
            end
        else
            push!(new_constraints_list, constraint)
        end
    end
    # 处理客体
    for constraint in another_top.constraints
        if length(over_a1) != 0 && length(over_a2) != 0
            if current_first                                                # 排除客体A2-A2
                if !(constraint.ai.nr in over_a2) && !(constraint.aj.nr in over_a2)
                    push!(new_constraints_list, constraint)
                end
            else                                                            # 排除客体A2-A2, 客体A2-B2替换为A1-B2
                if !(constraint.ai.nr in over_a2) || !(constraint.aj.nr in over_a2)
                    if constraint.ai.nr in over_a2
                        constraint.ai = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== constraint.ai.nr)])
                    end
                    if constraint.aj.nr in over_a2
                        constraint.aj = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== constraint.aj.nr)])
                    end
                    push!(new_constraints_list, constraint)
                end
            end
        else
            push!(new_constraints_list, constraint)
        end
    end
    top.constraints = new_constraints_list

    # Combine angles part
    println("Appending angles of $(another_top.moleculetype) to $(top.moleculetype) ...")
    new_angles_list = []
    # 处理主体
    for angle in top.angles
        if length(over_b1) != 0 && length(over_b2) != 0
            if current_first                                # 排除主体B1-B1, 主体A1-B1替换成A1-B2
                if !(angle.ai.nr in over_b1) || !(angle.aj.nr in over_b1) || !(angle.ak.nr in over_b1)
                    if angle.ai.nr in over_b1
                        angle.ai = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== angle.ai.nr)])
                    end
                    if angle.aj.nr in over_b1
                        angle.aj = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== angle.aj.nr)])
                    end
                    if angle.ak.nr in over_b1
                        angle.ak = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== angle.ak.nr)])
                    end
                    push!(new_angles_list, angle)
                end
            else                                            # 排除主体A1-B1
                if !(angle.ai.nr in over_b1) && !(angle.aj.nr in over_b1) && !(angle.ak.nr in over_b1)
                    push!(new_angles_list, angle)
                end
            end
        else
            push!(new_angles_list, angle)
        end
    end
    # 处理客体
    for angle in another_top.angles
        if length(over_a1) != 0 && length(over_a2) != 0
            if current_first                                                # 排除客体A2-A2
                if !(angle.ai.nr in over_a2) && !(angle.aj.nr in over_a2) && !(angle.ak.nr in over_a2)
                    push!(new_angles_list, angle)
                end
            else                                                            # 排除客体A2-A2, 客体A2-B2替换为A1-B2
                if !(angle.ai.nr in over_a2) || !(angle.aj.nr in over_a2) || !(angle.ak.nr in over_a2)
                    if angle.ai.nr in over_a2
                        angle.ai = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== angle.ai.nr)])
                    end
                    if angle.aj.nr in over_a2
                        angle.aj = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== angle.aj.nr)])
                    end
                    if angle.ak.nr in over_a2
                        angle.ak = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== angle.ak.nr)])
                    end
                    push!(new_angles_list, angle)
                end
            end
        else
            push!(new_angles_list, angle)
        end
    end
    top.angles = new_angles_list

    # Combine dihedrals part
    println("Appending dihedrals of $(another_top.moleculetype) to $(top.moleculetype) ...")
    new_dihedrals_list = []
    # 处理主体
    for dihedral in top.dihedrals
        if length(over_b1) != 0 && length(over_b2) != 0
            if current_first                                # 排除主体B1-B1, 主体A1-B1替换成A1-B2
                if !(dihedral.ai.nr in over_b1) || !(dihedral.aj.nr in over_b1) || !(dihedral.ak.nr in over_b1) || !(dihedral.al.nr in over_b1)
                    if dihedral.ai.nr in over_b1
                        dihedral.ai = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== dihedral.ai.nr)])
                    end
                    if dihedral.aj.nr in over_b1
                        dihedral.aj = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== dihedral.aj.nr)])
                    end
                    if dihedral.ak.nr in over_b1
                        dihedral.ak = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== dihedral.ak.nr)])
                    end
                    if dihedral.al.nr in over_b1
                        dihedral.al = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== dihedral.al.nr)])
                    end
                    push!(new_dihedrals_list, dihedral)
                end
            else                                            # 排除主体A1-B1
                if !(dihedral.ai.nr in over_b1) && !(dihedral.aj.nr in over_b1) && !(dihedral.ak.nr in over_b1) && !(dihedral.al.nr in over_b1)
                    push!(new_dihedrals_list, dihedral)
                end
            end
        else
            push!(new_dihedrals_list, dihedral)
        end
    end
    # 处理客体
    for dihedral in another_top.dihedrals
        if length(over_a1) != 0 && length(over_a2) != 0
            if current_first                                                # 排除客体A2-A2
                if !(dihedral.ai.nr in over_a2) && !(dihedral.aj.nr in over_a2) && !(dihedral.ak.nr in over_a2) && !(dihedral.al.nr in over_a2)
                    push!(new_dihedrals_list, dihedral)
                end
            else                                                            # 排除客体A2-A2, 客体A2-B2替换为A1-B2
                if !(dihedral.ai.nr in over_a2) || !(dihedral.aj.nr in over_a2) || !(dihedral.ak.nr in over_a2) || !(dihedral.al.nr in over_a2)
                    if dihedral.ai.nr in over_a2
                        dihedral.ai = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== dihedral.ai.nr)])
                    end
                    if dihedral.aj.nr in over_a2
                        dihedral.aj = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== dihedral.aj.nr)])
                    end
                    if dihedral.ak.nr in over_a2
                        dihedral.ak = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== dihedral.ak.nr)])
                    end
                    if dihedral.al.nr in over_a2
                        dihedral.al = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== dihedral.al.nr)])
                    end
                    push!(new_dihedrals_list, dihedral)
                end
            end
        else
            push!(new_dihedrals_list, dihedral)
        end
    end
    top.dihedrals = new_dihedrals_list

    # Combine exclusions part
    println("Appending exclusions of $(another_top.moleculetype) to $(top.moleculetype) ...")
    new_exclusions_list = []
    # 处理主体
    for exclusion in top.exclusions
        if length(over_b1) != 0 && length(over_b2) != 0
            if current_first                                # 排除主体B1-B1, 主体A1-B1替换成A1-B2
                if !(exclusion.ai.nr in over_b1) || !(exclusion.aj.nr in over_b1)
                    if exclusion.ai.nr in over_b1
                        exclusion.ai = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== exclusion.ai.nr)])
                    end
                    if exclusion.aj.nr in over_b1
                        exclusion.aj = get_atom_from_nr(another_top.atoms, over_b2[findfirst(over_b1 .== exclusion.aj.nr)])
                    end
                    push!(new_exclusions_list, exclusion)
                end
            else                                            # 排除主体A1-B1
                if !(exclusion.ai.nr in over_b1) && !(exclusion.aj.nr in over_b1)
                    push!(new_exclusions_list, exclusion)
                end
            end
        else
            push!(new_exclusions_list, exclusion)
        end
    end
    # 处理客体
    for exclusion in another_top.exclusions
        if length(over_a1) != 0 && length(over_a2) != 0
            if current_first                                                # 排除客体A2-A2
                if !(exclusion.ai.nr in over_a2) && !(exclusion.aj.nr in over_a2)
                    push!(new_exclusions_list, exclusion)
                end
            else                                                            # 排除客体A2-A2, 客体A2-B2替换为A1-B2
                if !(exclusion.ai.nr in over_a2) || !(exclusion.aj.nr in over_a2)
                    if exclusion.ai.nr in over_a2
                        exclusion.ai = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== exclusion.ai.nr)])
                    end
                    if exclusion.aj.nr in over_a2
                        exclusion.aj = get_atom_from_nr(top.atoms, over_a1[findfirst(over_a2 .== exclusion.aj.nr)])
                    end
                    push!(new_exclusions_list, exclusion)
                end
            end
        else
            push!(new_exclusions_list, exclusion)
        end
    end
    top.exclusions = new_exclusions_list

    # Combine atoms part
    shift_residues!(another_top, top.atoms[end].resnr)  # 客体残基编号平移
    println("Appending atoms of $(another_top.moleculetype) to $(top.moleculetype) ...")
    new_atoms_list = []
    for atom in top.atoms
        if !(atom.nr in over_b1)
            push!(new_atoms_list, atom)
        end
    end
    for atom in another_top.atoms
        if !(atom.nr in over_a2)
            push!(new_atoms_list, atom)
        end
    end
    top.atoms = new_atoms_list

    resort!(top.atoms)

    finish_hint("Finished appending $(another_top.moleculetype) to $(top.moleculetype).\n")
    title_hint("Write to itp now (option 0) or do other operations.\n")
end

end # module
